out <- readRDS("depends/out.rds")

mcmc <- out$mcmc

bayesplot::color_scheme_set("viridisA")
ggplot2::theme_set(theme_minimal())

#' Rhat
#' Looking for values less than 1.1 here
rhats <- bayesplot::rhat(mcmc)
bayesplot::mcmc_rhat(rhats)

(big_rhats <- rhats[rhats > 1.1])
##          beta_rho[1]          beta_rho[2]        beta_alpha[1]        beta_alpha[2]      logit_phi_rho_x      log_sigma_rho_x     logit_phi_rho_xs     log_sigma_rho_xs      logit_phi_rho_a      log_sigma_rho_a     logit_phi_rho_as     log_sigma_rho_as           u_rho_x[4] 
##             1.436462             2.050365             1.273354             1.462126             1.248576             1.150130             1.494549             1.362138             1.126501             1.151890             1.256854             1.325939             1.171248 
##          us_rho_x[4]          us_rho_x[6]          u_rho_xs[1]          u_rho_xs[2]          u_rho_xs[4]         u_rho_xs[11]         u_rho_xs[12]         u_rho_xs[22]         u_rho_xs[27]         u_rho_xs[30]         us_rho_xs[1]         us_rho_xs[2]         us_rho_xs[3] 
##             1.139724             1.159936             1.108818             1.104644             1.118410             1.113635             1.110012             1.174033             1.201293             1.139680             1.127415             1.138374             1.174318 
##         us_rho_xs[4]         us_rho_xs[5]         us_rho_xs[6]         us_rho_xs[8]         us_rho_xs[9]        us_rho_xs[10]        us_rho_xs[11]        us_rho_xs[12]        us_rho_xs[13]        us_rho_xs[14]        us_rho_xs[15]        us_rho_xs[22]        us_rho_xs[23] 
##             1.133138             1.156225             1.183878             1.133482             1.202853             1.187893             1.199672             1.198588             1.185517             1.125922             1.194364             1.185790             1.110124 
##        us_rho_xs[24]        us_rho_xs[25]        us_rho_xs[26]        us_rho_xs[27]        us_rho_xs[28]        us_rho_xs[29]        us_rho_xs[30]        us_rho_xs[31]        us_rho_xs[32]           u_rho_a[1]           u_rho_a[2]           u_rho_a[3]           u_rho_a[4] 
##             1.238823             1.134849             1.259984             1.296719             1.176014             1.247366             1.324837             1.255172             1.188287             1.429221             1.429893             1.435065             1.430681 
##           u_rho_a[5]           u_rho_a[6]           u_rho_a[7]           u_rho_a[8]           u_rho_a[9]          u_rho_a[10]          u_rho_as[1]          u_rho_as[2]          u_rho_as[3]          u_rho_as[4]          u_rho_as[5]          u_rho_as[6]          u_rho_as[7] 
##             1.432134             1.432472             1.432241             1.436418             1.434241             1.434673             2.017298             2.033510             2.041822             2.038281             2.054323             2.044709             2.034816 
##          u_rho_as[8]          u_rho_as[9]         u_rho_as[10]    logit_phi_alpha_x    log_sigma_alpha_x   logit_phi_alpha_xs   log_sigma_alpha_xs         u_alpha_x[1]         u_alpha_x[2]         u_alpha_x[3]         u_alpha_x[4]         u_alpha_x[8]         u_alpha_x[9] 
##             2.029767             2.036252             1.996995             1.101527             1.410485             1.146635             1.268746             1.263213             1.166815             1.240560             1.177752             1.193106             1.171466 
##        u_alpha_x[10]        u_alpha_x[13]        u_alpha_x[14]        u_alpha_x[21]        u_alpha_x[22]        u_alpha_x[23]        u_alpha_x[24]        u_alpha_x[25]        u_alpha_x[26]        u_alpha_x[27]        u_alpha_x[28]        u_alpha_x[30]        us_alpha_x[1] 
##             1.187447             1.273750             1.117343             1.432030             1.195067             1.190772             1.272183             1.363360             1.142784             1.178192             1.109855             1.130643             1.168862 
##        us_alpha_x[2]        us_alpha_x[3]        us_alpha_x[4]        us_alpha_x[5]        us_alpha_x[6]        us_alpha_x[8]       us_alpha_x[12]       us_alpha_x[15]       us_alpha_x[21]       us_alpha_x[22]       us_alpha_x[24]       us_alpha_x[25]       us_alpha_x[29] 
##             1.133704             1.197043             1.264392             1.112401             1.198646             1.100455             1.118491             1.118520             1.189470             1.245977             1.164936             1.132086             1.162871 
##       us_alpha_x[30]        u_alpha_xs[1]        u_alpha_xs[8]       u_alpha_xs[11]       u_alpha_xs[13]       u_alpha_xs[16]       u_alpha_xs[17]       u_alpha_xs[18]       u_alpha_xs[19]       u_alpha_xs[21]       u_alpha_xs[22]       u_alpha_xs[23]       u_alpha_xs[24] 
##             1.136288             1.115010             1.139084             1.103883             1.162435             1.145249             1.133684             1.299266             1.209079             1.405892             1.201692             1.144698             1.275732 
##       u_alpha_xs[25]       u_alpha_xs[28]       u_alpha_xs[29]       u_alpha_xs[32]       us_alpha_xs[1]       us_alpha_xs[2]       us_alpha_xs[3]       us_alpha_xs[4]       us_alpha_xs[5]       us_alpha_xs[6]       us_alpha_xs[8]       us_alpha_xs[9]      us_alpha_xs[10] 
##             1.257591             1.187026             1.390206             1.111991             1.158874             1.175500             1.192279             1.138272             1.181611             1.150340             1.217743             1.187316             1.192677 
##      us_alpha_xs[11]      us_alpha_xs[12]      us_alpha_xs[13]      us_alpha_xs[15]      us_alpha_xs[17]      us_alpha_xs[18]      us_alpha_xs[19]      us_alpha_xs[20]      us_alpha_xs[21]      us_alpha_xs[22]      us_alpha_xs[23]      us_alpha_xs[24]      us_alpha_xs[25] 
##             1.180764             1.122971             1.176139             1.124448             1.145104             1.138565             1.131682             1.104125             1.154424             1.203935             1.146553             1.133402             1.145399 
##      us_alpha_xs[26]         u_alpha_a[1]         u_alpha_a[2]         u_alpha_a[3]         u_alpha_a[4]         u_alpha_a[5]         u_alpha_a[6]         u_alpha_a[7]         u_alpha_a[8]         u_alpha_a[9]        u_alpha_a[10]        u_alpha_a[11]        u_alpha_a[12] 
##             1.115193             1.237430             1.260247             1.265380             1.224632             1.234161             1.248809             1.257960             1.259736             1.262103             1.253853             1.255800             1.234374 
##        u_alpha_a[13]        u_alpha_as[1]        u_alpha_as[2]        u_alpha_as[3]        u_alpha_as[4]        u_alpha_as[5]        u_alpha_as[6]        u_alpha_as[7]        u_alpha_as[8]        u_alpha_as[9]       u_alpha_as[10]       u_alpha_xa[13]       u_alpha_xa[21] 
##             1.218461             1.405071             1.431038             1.426739             1.435003             1.462654             1.469958             1.488702             1.483664             1.495276             1.430269             1.107958             1.335780 
##       u_alpha_xa[24]       u_alpha_xa[25]       u_alpha_xa[29]   log_sigma_lambda_x       ui_lambda_x[1]      ui_lambda_x[10]      ui_lambda_x[11]      ui_lambda_x[13]      ui_lambda_x[14]      ui_lambda_x[22]      ui_lambda_x[23]      ui_lambda_x[32] log_sigma_ancalpha_x 
##             1.259357             1.108961             1.204265             2.309844             1.212403             1.427155             1.258550             1.187641             1.158190             1.139450             1.131882             1.146537             1.470755 
##     ui_anc_rho_x[15]     ui_anc_rho_x[25]     ui_anc_rho_x[29]   ui_anc_alpha_x[16]   ui_anc_alpha_x[18]   ui_anc_alpha_x[19]   ui_anc_alpha_x[24]   ui_anc_alpha_x[25]   ui_anc_alpha_x[29]   ui_anc_alpha_x[32]      log_or_gamma[3]      log_or_gamma[5]      log_or_gamma[8] 
##             1.101091             1.172969             1.114979             1.195326             1.303179             1.144301             1.125907             1.151146             1.269189             1.129052             1.124149             1.113333             1.101212 
##      log_or_gamma[9]     log_or_gamma[10]     log_or_gamma[12]     log_or_gamma[18]     log_or_gamma[20]     log_or_gamma[22]     log_or_gamma[31]                 lp__ 
##             1.123776             1.127447             1.141006             1.102396             1.114067             1.128760             1.113397             1.567209
length(big_rhats) / length(rhats)
## [1] 0.4126016
#' ESS ratio
#' Worry about values less than 0.1 here... which they all are :clown:
ratios <- bayesplot::neff_ratio(mcmc)
bayesplot::mcmc_neff(ratios)

#' Autocorrelation
#' High values of autocorrelation in the chains
bayesplot::mcmc_acf(mcmc, pars = vars(starts_with("beta")))

#' Univariate traceplots

#' Assess these visually
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("beta")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("logit")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_sigma")))

#' Prevalence model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_a[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_as[")))

#' ART model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_a[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_as[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xa[")))

#' Other
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_lambda_x[")))

#' ANC model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_or_gamma["))) #' N.B. these are from the ANC attendance model

#' Pairs plots
#' Prior suspicion (from Jeff, Tim, Rachel) that the ART attendance model is unidentifiable
#' Let's have a look at the pairs plot for neighbouring districts and the log_or_gamma parameter
nb <- area_merged %>%
  filter(area_level == max(area_level)) %>%
  bsae::sf_to_nb()

neighbours_pairs_plot <- function(par, i) {
  neighbour_pars <- paste0(par, "[", c(i, nb[[i]]), "]")
  bayesplot::mcmc_pairs(mcmc, pars = neighbour_pars, diag_fun = "hist", off_diag_fun = "hex")
}

area_merged %>%
  filter(area_level == max(area_level)) %>%
  print(n = Inf)
## Simple feature collection with 32 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 32.67152 ymin: -17.12721 xmax: 35.91505 ymax: -9.364003
## Geodetic CRS:  WGS 84
## # A tibble: 32 × 12
##    area_id       area_name     area_level parent_area_id spectrum_region_code area_sort_order center_x center_y area_level_label display naomi_level                                                                                geometry
##  * <chr>         <chr>              <int> <chr>                         <dbl>           <int>    <dbl>    <dbl> <chr>            <lgl>   <lgl>                                                                            <MULTIPOLYGON [°]>
##  1 MWI_4_1_demo  Chitipa                4 MWI_3_1_demo                      0              38     33.4    -9.98 District + Metro TRUE    TRUE        (((33.63988 -9.603465, 33.64077 -9.611695, 33.64722 -9.620278, 33.65034 -9.630591, 3...
##  2 MWI_4_2_demo  Karonga                4 MWI_3_2_demo                      0              39     33.9   -10.1  District + Metro TRUE    TRUE        (((34.17391 -10.58124, 34.16676 -10.57488, 34.15757 -10.57275, 34.15369 -10.56901, 3...
##  3 MWI_4_3_demo  Rumphi                 4 MWI_3_3_demo                      0              40     33.8   -10.8  District + Metro TRUE    TRUE        (((33.38671 -11.14385, 33.37861 -11.14527, 33.37376 -11.13834, 33.37219 -11.12912, 3...
##  4 MWI_4_4_demo  Mzuzu City             4 MWI_3_4_demo                      0              41     34.0   -11.4  District + Metro TRUE    TRUE        (((33.96686 -11.49152, 33.96384 -11.48679, 33.95983 -11.483, 33.95649 -11.47607, 33....
##  5 MWI_4_5_demo  Nkhata Bay             4 MWI_3_5_demo                      0              42     34.0   -11.7  District + Metro TRUE    TRUE        (((34.01984 -12.23471, 34.01532 -12.23359, 34.00829 -12.2278, 34.0042 -12.22929, 33....
##  6 MWI_4_6_demo  Mzimba                 4 MWI_3_4_demo                      0              43     33.6   -11.8  District + Metro TRUE    TRUE        (((33.53559 -12.37776, 33.54323 -12.37243, 33.54602 -12.36943, 33.54572 -12.36572, 3...
##  7 MWI_4_7_demo  Likoma                 4 MWI_3_6_demo                      0              44     34.7   -12.1  District + Metro TRUE    TRUE        (((34.73153 -12.03097, 34.73347 -12.03236, 34.73292 -12.03681, 34.73486 -12.03875, 3...
##  8 MWI_4_8_demo  Nkhotakota             4 MWI_3_7_demo                      0              45     34.1   -12.8  District + Metro TRUE    TRUE        (((34.12471 -12.3983, 34.13183 -12.39831, 34.14007 -12.3996, 34.14636 -12.40232, 34....
##  9 MWI_4_9_demo  Kasungu                4 MWI_3_8_demo                      0              46     33.5   -13.0  District + Metro TRUE    TRUE        (((33.32294 -13.61568, 33.31926 -13.61632, 33.31408 -13.6193, 33.3055 -13.61247, 33....
## 10 MWI_4_10_demo Ntchisi                4 MWI_3_9_demo                      0              47     33.9   -13.3  District + Metro TRUE    TRUE        (((34.1272 -13.50292, 34.12042 -13.50349, 34.11702 -13.50761, 34.1156 -13.51225, 34....
## 11 MWI_4_11_demo Dowa                   4 MWI_3_10_demo                     0              48     33.7   -13.5  District + Metro TRUE    TRUE        (((33.83662 -13.7742, 33.83127 -13.77034, 33.82476 -13.75987, 33.82003 -13.75745, 33...
## 12 MWI_4_12_demo Salima                 4 MWI_3_11_demo                     0              49     34.4   -13.8  District + Metro TRUE    TRUE        (((34.55063 -14.10458, 34.54295 -14.10673, 34.53391 -14.1058, 34.52168 -14.10958, 34...
## 13 MWI_4_13_demo Mchinji                4 MWI_3_12_demo                     0              50     33.0   -13.8  District + Metro TRUE    TRUE        (((32.92279 -13.37107, 32.92497 -13.36179, 32.92436 -13.3568, 32.94433 -13.35478, 32...
## 14 MWI_4_14_demo Lilongwe City          4 MWI_3_13_demo                     0              51     33.8   -13.9  District + Metro TRUE    TRUE        (((33.75722 -13.75277, 33.76022 -13.75845, 33.77182 -13.76671, 33.77989 -13.76876, 3...
## 15 MWI_4_15_demo Lilongwe               4 MWI_3_13_demo                     0              52     33.5   -14.0  District + Metro TRUE    TRUE        (((33.70239 -14.54722, 33.71319 -14.57181, 33.70703 -14.57638, 33.70465 -14.58096, 3...
## 16 MWI_4_16_demo Dedza                  4 MWI_3_14_demo                     0              53     34.3   -14.2  District + Metro TRUE    TRUE        (((34.40013 -14.412, 34.39294 -14.40303, 34.39324 -14.3956, 34.38913 -14.39334, 34.3...
## 17 MWI_4_17_demo Ntcheu                 4 MWI_3_15_demo                     0              54     34.8   -14.8  District + Metro TRUE    TRUE        (((34.80996 -15.32603, 34.80709 -15.32502, 34.8033 -15.3176, 34.80719 -15.31364, 34....
## 18 MWI_4_18_demo Mangochi               4 MWI_3_16_demo                     0              55     35.3   -14.1  District + Metro TRUE    TRUE        (((34.97112 -14.77456, 34.96737 -14.7704, 34.95998 -14.75692, 34.95798 -14.75167, 34...
## 19 MWI_4_19_demo Machinga               4 MWI_3_17_demo                     0              56     35.6   -14.9  District + Metro TRUE    TRUE        (((35.16126 -15.20686, 35.16416 -15.19912, 35.16824 -15.19152, 35.1692 -15.18516, 35...
## 20 MWI_4_20_demo Balaka                 4 MWI_3_18_demo                     0              57     35.1   -15.0  District + Metro TRUE    TRUE        (((35.03931 -15.3229, 35.00434 -15.32357, 34.98965 -15.3241, 34.91631 -15.32533, 34....
## 21 MWI_4_21_demo Zomba City             4 MWI_3_19_demo                     0              58     35.3   -15.4  District + Metro TRUE    TRUE        (((35.34807 -15.41108, 35.34317 -15.4197, 35.3424 -15.42337, 35.33797 -15.42816, 35....
## 22 MWI_4_22_demo Zomba                  4 MWI_3_19_demo                     0              59     35.4   -15.4  District + Metro TRUE    TRUE        (((35.84263 -15.45047, 35.83955 -15.44715, 35.83957 -15.44162, 35.83749 -15.43646, 3...
## 23 MWI_4_23_demo Phalombe               4 MWI_3_20_demo                     0              60     35.7   -15.7  District + Metro TRUE    TRUE        (((35.80985 -15.93286, 35.80041 -15.93066, 35.79598 -15.93278, 35.78529 -15.92961, 3...
## 24 MWI_4_24_demo Mulanje                4 MWI_3_21_demo                     0              61     35.5   -15.9  District + Metro TRUE    TRUE        (((35.30769 -16.21557, 35.30645 -16.21306, 35.29699 -16.20214, 35.28732 -16.19913, 3...
## 25 MWI_4_25_demo Neno                   4 MWI_3_22_demo                     0              62     34.7   -15.5  District + Metro TRUE    TRUE        (((34.73322 -15.81642, 34.72392 -15.81112, 34.71292 -15.81164, 34.71111 -15.80616, 3...
## 26 MWI_4_26_demo Blantyre               4 MWI_3_23_demo                     0              63     34.9   -15.7  District + Metro TRUE    TRUE        (((34.90176 -16.01093, 34.8995 -16.01271, 34.89635 -16.00977, 34.8893 -16.01434, 34....
## 27 MWI_4_27_demo Mwanza                 4 MWI_3_24_demo                     0              64     34.5   -15.6  District + Metro TRUE    TRUE        (((34.73546 -15.81671, 34.73438 -15.81708, 34.52702 -15.81355, 34.52308 -15.80903, 3...
## 28 MWI_4_28_demo Chiradzulu             4 MWI_3_25_demo                     0              65     35.2   -15.8  District + Metro TRUE    TRUE        (((35.31332 -15.98932, 35.30383 -16.00106, 35.30011 -16.00326, 35.2948 -16.00073, 35...
## 29 MWI_4_29_demo Blantyre City          4 MWI_3_23_demo                     0              66     35.1   -15.8  District + Metro TRUE    TRUE        (((35.10285 -15.85812, 35.1012 -15.86633, 35.09855 -15.87003, 35.09046 -15.87081, 35...
## 30 MWI_4_30_demo Thyolo                 4 MWI_3_26_demo                     0              67     35.1   -16.1  District + Metro TRUE    TRUE        (((35.26228 -16.40585, 35.25777 -16.4048, 35.25433 -16.40626, 35.24886 -16.4003, 35....
## 31 MWI_4_31_demo Chikwawa               4 MWI_3_27_demo                     0              68     34.7   -16.2  District + Metro TRUE    TRUE        (((34.9105 -16.69166, 34.90701 -16.68995, 34.88869 -16.68457, 34.8867 -16.6869, 34.8...
## 32 MWI_4_32_demo Nsanje                 4 MWI_3_28_demo                     0              69     35.1   -16.7  District + Metro TRUE    TRUE        (((34.9105 -16.69166, 34.91037 -16.688, 34.90659 -16.68455, 34.90675 -16.67798, 34.9...
neighbours_pairs_plot("log_or_gamma", 5) #' Nkhata Bay and neighbours

neighbours_pairs_plot("log_or_gamma", 26) #' Blantyre and neighbours

#' NUTS specific assessment
np <- bayesplot::nuts_params(mcmc)

#' Number of divergent transitions
#' In this instance it's zero, so no need to do further divergent transition investigation
np %>%
  filter(Parameter == "divergent__") %>%
  summarise(n_divergent = sum(Value))
bayesplot::mcmc_nuts_divergence(np, bayesplot::log_posterior(mcmc))

#' Energy plots
#' Ideally these two histograms would be the same (Betancourt 2017)
#' The histograms are quite different, in a way suggesting the chains may not be
#' fully exploring the tails of the target distribution
bayesplot::mcmc_nuts_energy(np)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.